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Abstract 

We study the numerical solution of the time-dependent Gross-Pitaevskii 
equation (GPE) describing a Bose-Einstein condensate (BEC) at zero or very 
low temperature. In preparation for the numerics we scale the 3d Gross- 
Pitaevskii equation and obtain a four-parameter model. Identifying 'extreme 
parameter regimes', the model is accessible to analytical perturbation theory, 
which justifies formal procedures well known in the physical literature: re- 
duction to 2d and Id GPEs, approximation of ground state solutions of the 
GPE and geometrical optics approximations. Then we use a time-splitting 
spectral method to discretize the time-dependent GPE. Again, perturbation 
theory is used to understand the discretization scheme and to choose the spa- 
tial/temporal grid in dependence of the perturbation parameter. Extensive 
numerical examples in Id, 2d and 3d for weak/strong interactions, defocus- 
ing/focusing nonlinearity, and zero/nonzero initial phase data are presented 
to demonstrate the power of the numerical method and to discuss the physics 
of Bose-Einstein condensation. 

*Email address: bao@cz3.nus.edu.sg. Fax: 65-67746756 
tEmail address: d.jakschl@physics.ox.ac.uk. 

■^Email address: peter. markowich@univie.ac. at, http://mailbox.univie.ac.at/peter.markowich. 



1 



Key Words: Bose-Einstein condensation (BEC), Gross-Pitaevskii equation, time- 
splitting spectral method, approximate ground state solution, defocusing/focusing 
nonlinearity. 

1 Introduction 

Recent experimental advances in achieving and observing Bose-Einstein condensa- 
tion (BEC) in trapped neutral atomic vapors [5, 22, 12] have spurred great excite- 
ment in the atomic physics community and renewed the interest in studying the 
collective dynamics of macroscopic ensembles of atoms occupying the same one- 
particle quantum state [46, 21, 33]. The condensate typically consists of a few 
thousands to millions of atoms which are confined by a trap potential. In fact, 
beside the effects of the internal interactions between the atoms, the macroscopic 
behavior of BEC matter is highly sensitive to the shape of this external trapping 
potential. Theoretical predictions of the properties of a BEC like the density profile 
[11], collective excitations [25] and the formation of vortices [50] can now be com- 
pared with experimental data [5, 37, 45]. Needless to say that this dramatic progress 
on the experimental front has stimulated a wave of activity on both the theoretical 
and the numerical front. 

The properties of a BEC at temperatures T much smaller than the critical con- 
densation temperature T c [39] are usually well modeled by a nonlinear Schrodinger 
equation (NLSE) for the macroscopic wave function [39, 33] known as the Gross- 
Pitaevskii equation (GPE) [34, 48], which incorporates the trap potential as well as 
the interactions among the atoms. The effect of the interactions is described by a 
mean field which leads to a nonlinear term in the GPE. The cases of repulsive and 
attractive interactions - which can both be realized in the experiment - correspond 
to defocusing and focusing nonlinearities in the GPE, respectively. Note that equa- 
tions very similar to the GPE also appear in nonlinear optics where an index of 
refraction, which depends on the light intensity, leads to a nonlinear term like the 
one encountered in the GPE. 

There has been a series of recent studies which deals with the numerical solution 
of the time-independent GPE for the ground state and the time-dependent GPE for 
finding the dynamics of a BEC. Bao et al. [9] presented a general method to compute 
the ground state solution via directly minimizing the energy functional and used it to 
compute the ground state of the GPE in different cases. Edwards et al. introduced a 
Runge-Kutta type method and employed it to solve the spherically symmetric time- 
independent GPE [26]. Adhikari [1, 2] used this approach to obtain the ground state 
solution of the GPE in 2d with radial symmetry. Other approaches include a finite 
difference method proposed by Chiofalo et al. [18] and a simple analytical method 
proposed by Dodd [23] . For the numerical solution of the time-dependent GPE only 
few methods are available, a particle-inspired scheme proposed by Cerimele et al. in 
[16] and a finite difference method used by Ruprecht et al. [51]. 

In this paper, we take the 3d Gross-Pitaevskii equation, make it dimensionless to 
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obtain a four-parameter model, use (singular) perturbation theory to discuss semi- 
classical asymptotics, to approximately reduce it to a 2d GPE and a Id GPE in 
certain limits, and discuss the approximate ground state solution of the GPE in two 
extreme regimes: (very) weak interactions and strong repulsive interactions, again 
using perturbation methods. Numerical computations for similar physical set ups 
can to a large extent be found in the physical literatures cf. [11], however they are 
included here since perturbation theory gives a systematic way to obtain rigorously 
(validate) these approximations and since they are used as important preparatory 
steps for the numerical simulations (approximate ground states usually serve as ini- 
tial data). Then we use the time-splitting spectral method, which was studied in 
Bao et al. [7, 8] for the Schrodinger equation in the semiclassical regime, to discretize 
the time-dependent GPE. The merit of the numerical method is that it is explicit, 
unconditionally stable, time reversible, time-transverse invariant, and conserves the 
position density. In fact, the spectral method has shown great success in solving 
problems arising from many areas of science [32, 13] due to its spatially spectral ac- 
curacy. The split-step procedure was presented for differential equations in [54] and 
applied to Schrodinger equations [35, 55, 28] and the KDV equation [56], as well as 
used with an iterative procedure for optical fibers [4]. In this paper we perform grid 
control using singular perturbation theory and present extensive numerical examples 
in Id, 2d and 3d for weak/strong interactions, defocusing/focusing nonlinearity, and 
zero/nonzero initial phase data. 

The paper is organized as follows. In section 2 we start out with the 3d GPE, 
scale it to get a four-parameter model, show how to reduce it to lower dimensions and 
give the approximate ground state solution in the two mentioned extreme regimes of 
(very) weak interactions and strong repulsive interactions and discuss semiclassical 
asymptotics. In section 3 we present the time-splitting spectral method for the GPE. 
In section 4 numerical tests of the GPE for different cases including weak/strong 
interactions, defocusing/focusing nonlinearity, and zero/nonzero initial phase data 
are presented. In section 5 a summary is given. 



2 Gross-Pitaevskii equation 

At temperatures T much smaller than the critical temperature T c [39] , a BEC is well 
described by the macroscopic wave function ip — ^(x, t) whose evolution is governed 
by a self-consistent, mean field nonlinear Schrodinger equation (NLSE) known as the 
Gross-Pitaevskii equation [34, 48]. If a harmonic trap potential is considered, the 
equation becomes 

lh d J^A = ~V>(x,t)+| (^+,y + W2 V) ^x,t)+NU \^,t)\^(x,t), 

(2.1) 

where x = (x, y, z) T is the spatial coordinate vector, m is the atomic mass, h is 
the Planck constant, N is the number of atoms in the condensate, and u> x , uj v and 
uj z are the trap frequencies in x, y and z-direction, respectively. For the following 
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we assume (w.r.o.g.) oo x < oo y < uj z . When the trap potential is 

isotropic. U$ describes the interaction between atoms in the condensate and has the 
form 

U„ = (2.2) 
m 

where a is the s-wave scattering length (positive for repulsive interaction and neg- 
ative for attractive interaction). It is necessary to ensure that the wave function is 
properly normalized. Specifically, we require 

/ HUM 2 dx= 1. (2.3) 

JR 3 

A typical set of parameters used in current experiments with 87 Rb is given by 

m = 1 .44 xlCT 25 [£;#], u x = u y = to z = 20n[rad/s], a = 5.1xHT 9 [m], N : 10 2 ~ 10 7 

(2.4) 

and the Planck constant has the value 

% = 1.05 x 1(T 34 [Js]. 

2.1 Dimensionless GPE 

In order to scale the equation (2.1) under the normalization (2.3), we introduce 

i = u x t, x=— , ^(x,t) =x 3/ V(x,t), (2.5) 

x s 

where x s is the characteristic 'length' of the condensate. Plugging (2.5) into (2.1), 
multiplying by — rm, and then removing all ~, we obtain the following dimen- 

sionless Gross-Pitaevskii equation under the normalization (2.3) in three spatial 
dimensions 

l£ d J^A = _^v 2 ^(x,t) + V(xMx,t) + fe 5 / 2 |^(x,t)| 2 ^(x,t), (2.6) 



where 



nx) = i(x 2 + 7 y+7,v), 



n 


-m 


uj x mx 2 s 




U N 


4iraN 




a 



ly = i lz = 

Ul x 



5 = -^t — = , oo = 



with a the length of the harmonic oscillator ground state (in x-direction). The 
coefficient of the nonlinearity of (2.6) (interaction strength parameter) can also be 
expressed as 

K . = £ £ 5/2 = ^aN /ao\ 5 _ 1 8iraN a 4 = sgn(a) a 2 , a 2 , = sgn(a) /ao ao\ 2 , _v 
a \x s J 2 xj x 2 s 2 x\x\ 2 \x h x s J 




where Xh is the healing length [11] with 

Tr\n\N\ 2 

(2.8) 

If we plug the typical set of parameter values (2.4) into the above parameters, we 
find 

a « 0.3407 x l(T 5 [m], 5 w 0.01881iV : 1.881 - 188100. 

Remark 2.1 If one chooses x s = a in (2.5), then e = 1 in (2.6), k — 5 in 
(2.7), and the equation (2.6) takes the form often appearing in the physical litera- 
ture. This choice for x s is appropriate in the weak interaction regime characterized 
by 47r|a|iV <C ao and in the moderate interaction regime where 4ir\a\N ^ a . In 
the strong interaction regime 47r|a|iV ^> ao a different choice is more appropriate, 

namely x s = (47r|a|iVao) 1//5 , which gives \k\ = 1 and e = ( ^ffpv ) ^ ^ Other 
choices for x s , based on approximating the actual condensate length scale, shall be 
discussed in Section 2.2. Note that the choice of x s determines the observation 
scale of the condensate and decides: (i) which phenomena are 'visible' by asymp- 
totic analysis, (ii) which phenomena can be resolved by discretization on specified 
spatial/temporal grids. 



Thus, there are two extreme regimes: one is when e = 0(1) ao = 0(x s )) and 
k = fe 5 / 2 = o(l) 47r|a|iV <C ao), then equation (2.6) describes a weakly inter- 
acting condensate. The other is when e = o(l) x s ^> a ) and k = 5e 5 ^ 2 = 0(1) 
(implying 47r|a|iV ^> a ) (or e = 1 and k = 5e 5 ^ 2 = 5 with \S\ ^> 1 by the rescaling 
x — > £ 1 / 2 x, ip — > ip/e 3 ^), then (2.6) corresponds to a strongly interacting condensate 
(Thomas- Fermi regime [11]) or, equivalently, to the semiclassical regime. We recall 
that the equation (2.6) is regularly perturbed in the case of weak interactions and 
singularly perturbed in the semiclassical regime. Analytical techniques of asymp- 
totic analysis are available in both cases providing structural information on the 
solutions of (2.6) and on numerical discretization schemes (spatial/temporal grid 
control, control of the computational domain, error estimates in linearized cases). 

2.2 Approximate ground state solution in 3d 

To find a stationary state of (2.6), we write 

ip(x,t) = e-^' /£ 0(x), (2.9) 

where fj, is the chemical potential of the condensate. Inserting into (2.6) gives the 
following equation for 0(x) 

/i0(x) = ~ V 2 0(x) + V(x)0(x) + /#(x)| 2 0(x), x G K 3 , (2.10) 
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under the normalization condition 



(2.11) 



This is a nonlinear eigenvalue problem. The Bose-Einstein condensate ground-state 
wave function 9 (x) is found by solving this eigenvalue problem under the normal- 
ization condition (2.11) with the minimal chemical potential [i g . Usually, the ground 
state problem is formulated variationally. Define the energy functional 



It is easy to see that critical points of E are 'eigenfunctions' of the nonlinear Hamil- 
tonian. To compute the ground state <p g we solve 



In the case of a defocusing (stable) condensate the energy functional E(ip) is positive, 
coercive and weakly lower semicontinuous on the unit sphere in L 2 (R 3 ), thus the exis- 
tence of a minimum follows from standard theory. For understanding the uniqueness 
question note that E(acj) g ) = E(<p g ) for all a G C with \a\ = 1. Thus an additional 
constraint has to be introduced to show uniqueness, e.g. <f) g real valued and 9 (x) > 
for all x G IR 3 (see [41]). For focusing (unstable) 3-dimensional condensates the en- 
ergy functional E(ip) is not bounded from below on the unit sphere of L 2 (M 3 ). Thus, 
an absolute minimum of E(ip) does not exist on {ip G L 2 (R 3 ) | HV'HIz = 1}- The 
interpretation of critical points (local minimum, saddle points obtained by min-max- 
theory) as physically relevant ground states is by no means clear. 

Using (simple) perturbation methods, we present here the approximate ground 
state solution of (2.6) in the two extreme regimes of weak repulsive or attractive 
interactions and strong repulsive interactions (see [11] for a discussion in physical 
literature). 

These approximate ground state solutions are used in reducing the 3d GPE to a 
2d GPE and a Id GPE - see the next subsection for details - and as initial data for 
the numerical solution of the time- dependent GPE in section 4 (see the subsequent 
discussion). 

For a weakly interacting condensate, i.e. e = 0(1) and k = o(l), we drop the 
nonlinear term (i.e. the last term on the right hand side of (2.6)) and find the 
harmonic oscillator equation 



E(4)-= ^/ |V0| 2 rfx+ / U(x)|0| 2 rfx+^ / |0| 4 rfx. (2.12) 

2 JR3 JR3 2 JRi 



E{4> g ) 



min E((f>), 




(2.13) 



Mx) = -yV 2 0(x) + i (x 2 + 7 y + 72 V) 0(x). 



(2.14) 



The ground state solution of (2.14) is 




(2.15) 
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It can be viewed as an approximate ground state solution of (2.6) in the case of a 
weakly interacting condensate, with an 0(/t)-error in approximating the chemical 
potential and ground state. From (2.15), we can see that the diameter of the ground 
state solution (computed according to the formula (4.2)) in the weakly interacting 
condensate is 

x w g =0 (yi) = 0(1) (after the scaling (2.5)). 

Also we remark that the condensate widths in y and z-directions of the approximate 

ground state <j>™ are O = (Vv 7 ^) and (V^/V^) = (Vv 7 ^)' 

respectively. Clearly, this is important for the control of the computational domain. 

For a condensate with strong repulsive interactions, i.e. e = o(l), k = 0(1) and 
k > 0, we drop the diffusion term (i.e. the first term on the right hand side of (2.6)) 
corresponding to the Thomas- Fermi approximation [11]: 

^(x) = V(x)0(x) + /#(x)| 2 0(x), x G R 3 . (2.16) 

The ground state solution of (2.16) is the compactly supported function 



LS. 



e / 15^ 7z \ 2/5 _ 1 / 15^ 7 A 2 / 5 = | Mp* - V( X )) /«, V(x) < ft, 



2 V 47r / 2 V 4vr ; y o, otherwise. 

(2.17) 

This shows that the diameter of the ground state solution in the strongly interacting 
repulsive condensate is 

x s g= f^ g = 0((5 7s/ 7,) 1/5 )- = ( U*\a\N lylz \ /5 \ ao (aga . n after gcaling (2 5)) _ 

x s \ y flo J j x s 

Approximate widths of <j) s g in the y and ^-directions are O (W) V 77, 4/5 ) = O {ll^hT 

and O ((«7 J/ ) 1/5 /7 4/5 ) = O h A J b ) , respectively. 

This analysis suggests to choose the characteristic condensate length x s such 
that the (dimensionless) ground state width x™ is 0(1) after the scaling (2.5), i.e. 

x s = O(a ) for weak interaction (as in Remark 2.1) and x s = O ^( ^~ 7y7z) ^ a o^) 

for strong repulsive interaction (different from Remark 2.1 if ^ y ^> 1 or ^ z ^> 1). If 
we use the typical set of parameter values (2.4) in the above identity, we obtain 

e : 0.0078 ~ 1. 

Remark 2.2 The approximate ground state (p™ in the weak-interaction regime has 
finite energy, more precisely 

^W)=^ + l(^) 1/2 = ^ + °(«) ( 2 - 18 ) 
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(IS K > foT E = 0(1), lz = 0{l), ly = 0(1). 

Contrarily the energy of the Thomas- Fermi approximation is infinite 

E{cf> s g ) = +oo (2.19) 

due to the low regularity of <p s g at the free boundary V(x) = fi s g . More precisely, <f) s g is 
locally O 1 / 2 at the interface but not H\ oc (M. 3 ). This is a typical behavior for solutions 
of free boundary value problems, which indicates that <p s g does not approximate <p g 
to the full 0(e 2 ) - order, as indicated by formal consistency. An interface layer 
correction has to be constructed in order to improve the approximation quality. For 
a convergence proof of (j) s g — * (fi g (without convergence rate) we refer to [41]. 

It is of course tempting to use approximate ground states as initial data for the 
GPE when simulating Bose-Einstein condensation. In the weak interaction case 
this produces O(k) - errors in time dependent simulations on 0(1) time intervals. 
In the strong interaction case an initial wave function 4> s g produces time - dependent 
solutions with infinite energy (which usually generates breathing modes, cf. Example 
3 III in Sec. and the error in the wave function introduced by this is typically 
significantly larger than 0(e 2 ). 



2.3 Reduction to lower dimensions 

In two important cases, the 3d Gross-Pitaevskii equation (2.6) can approximately 
be reduced to a lower dimensional PDE. For a disk-shaped condensate with small 
height, i.e. 

lu x uj y , tu z > u x , j y ~ 1 , 7 Z > 1, (2.20) 

the 3d GPE (2.6) can be reduced to a 2d GPE with x = (x, y) T by assuming that the 
time evolution does not cause excitations along the z-axis since these have a large 
energy of approximately hu z compared to excitations along the x and y-axis with 
energies of about Tiuj x . To understand this, consider the total condensate energy 
E\ 



Em)\ = y / K3 iwwr^x + \ / r3 (x 2 + 7 y) m)\ 2 d* 

+ | j^mt)?^Jjm\^. (2.2D 

Multiplying (2.6) by ip t and integrating by parts show the energy conservation 

E[1>(t)\ = £[M Vt, (2.22) 

where ipj = ip(t = 0) is the initial function which may depend on all parameters e, 
7 y , 7 2 and k. Now assume that ipi satisfies 

-^^0, as 7,^oo. (2.23) 



7; 



z 



Take a sequence ;^ oo (and keep all other parameters Bxed). Since / |V-(i)| 2 <ix = 

1 we conclude from weak compactness that there is a positive measure n°(t) such 
that 

\i/j(t)\ 2 — 1 n°(t) weakly as 7 2 — > oo. 
Energy conservation implies 

/ z 2 |^(t) | 2 dx — > 0, as 7 2 — > oo 

Jr 3 

and thus we conclude concentration of the condensate in the plane z = 0: 

n°(x, y, z, t) = n° 2 (x, y, t)5(z), 

where n°(t) is a positive measure on R 2 . 

Now let fa = ipz(z) be a wave-function with 

/ \fo(z)\ 2 dz = l, 

JR 

depending on 7 2 such that 

\fo(z)\ 2 ^ 8(z), as 7 2 -*oo, (2.24) 
Denote by S the subspace 

s = {^ = ^,y)M^)\^eL 2 (R 2 )} 

and let 

IT : L 2 (R 3 ) -> S C L 2 (R 3 ) 

be the projection on 5: 

(n^)(x, y, z) = ip 3 (z) / $,(<r) tp(x, y, a) da. 
Now write the equation (2.6) in the form 

ieipt = Aip + 

where Aip stands for the linear part and J-{ijS) for the nonlinearity. Applying II to 
the GP-equation gives 

ie{Uip) t = ILAi/> + ILF(i/>) 

= iu(nv>) + iLF(ity) + n ((ru - An)^p + (ilf(^) - .F(ify))) .(2.25) 

The projection approximation of (2.6) is now obtained by dropping the commutator 
terms, it reads 

ie(Ua) t = TlA(Tla) + Oilier), (2.26) 
(IL7)(t = 0) = Ifyj, (2.27) 
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or explicitly, with 

(Ua)(x,y,z,t) =:^{x,y,t)ij} Z {z), (2.28) 

we find 

l£ ^f = -y v ^ 2 + 2 1 (* 2 + 7 ^ 2 + c ) ^ + ( fc5/2 /I ^ {z) dz ) (2 - 29) 

where 

/OO /"CO 
Z 2 |^W| 2 ^ + ^ 2 / 
-oo J — oc 

Since this GPE is time-transverse invariant, we can replace ip2 — ► ^ e~*fe and drop 
the constant C in the trap potential. The observables are not affected by this. 

The 'effective' GP-equation (2.29) is well known in the physical literature [40], 
where the projection method is often referred to as 'integrating out the z-coordinate'. 
However, an analysis of the limit process 72^00 has to be based on the derivation 
as presented above, in particular on studying the commutators IL4 — .All, ILF — FH. 
In the case of small interaction e = 0(1), k = o(l), a good choice for ^(z) is the 
ground state of the harmonic oscillator in z-dimension: 

^z) = g) V4 e -r.* 2 /&). (2.30) 

Note that |^ 3 (z)| 2 — 1 8{z) as 7^ — > 00 and that IIA = All such that the error in 
approximating Tlip(t) by Ua is determined by the commutator of the nonlinearity, 
which is O(k). 

For condensates with other than small interaction the choice of ^3 is much less 
obvious. Often one assumes that the condensate density along the z-axis well de- 
scribed by the (x, y)-trace of the ground state position density \<p g \ 2 

\ip(x,y,z,t)\ 2 « \^ 2 (x,y,t)\ 2 \(j) g (x 1 ,y 1 ,z)\ 2 dx^ (2.31) 

and (taking a pure-state-approximation) 

Mz) = (jf a \(f> g {x,y,z)\ 2 dxdy^ (2.32) 

A mathematical analysis of the limit process 7 2 — ■> 00 is currently under study. 

For a cigar-shaped condensate 

uj y > co x , u z > u x , 7 y > 1, 7 Z > 1, (2.33) 

the 3d GPE (2.6) can be reduced to a Id GPE by proceeding analogously. 

Then the 3d GPE (2.6), 2d GPE and Id GPE can then be written in a unified 

way 

l£ ^A = _^v 2 ^(x,t) + \/,(x)^(x,t) + ^ |#M)|VCM), x G R d , (2.34) 
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#3 



dz 



dz. 



where 



l„2 



X , 



' fvt?i/>2a(y,z) dydz, _ 

I K^ 2 +7y+7,v), 



The normalization condition for (2.34) is 

/ ^(x,t)| 2 dx=l. 



d = l, 
d = 2, 

d = 3. 
(2.35) 

(2.36) 



By using the approximate ground state of Section 2.2, we derive - after simple 
calculations - for a weakly interacting condensate 



n d := n d = < 



2tt£ 2vr ' 

° fc V27T' 



27T£ 



d=l, 
d = 2, 
d = 3, 



(2.37) 



and for a condensate with strong repulsive interactions 

8/5 



8/5 



i(sr^^7.)-/ 8 =^^(s) 
Kt) 1/5 ^=(^)^ 2 ( I fc) 1/5 ?, 



d= 1, 
d = 2, 
d = 3. 



(2.38) 



We call a rf-dimensional (d = 1 or 2) condensate 'very weakly interacting', if 
e = 0(1) and |/c d | = |/c™| <C 1, which implies I^I^/tJ «C 1 in 2d and |^|^/7 y 7z <C 1 
after reduction to Id. 



Remark 2.3 By using very elongated trapping potentials it is now possible to pro- 
duce weakly interacting condensates that are 'truly ' in the ID regime and the assump- 
tion that the condensate wave function factorizes is fulfilled with high accuracy. In 
cases where the interactions cannot be neglected low dimensional simulations [52] 
should be viewed as model calculations where the effective interaction strength in the 
reduced equation is estimated by Eq. (2.37) for a weakly interacting condensate or 
(2.38) for a condensate with strong repulsive interactions. More detailed studies go 
beyond the GPE to describe low dimensional EEC's [40]. 

2.4 Geometrical optics e — > 

We set 

^CM) = yjn&t) exp (-S(x,t) \ , 
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where n — \ip\ 2 and S is the phase of the wave-function. Inserting into the GP- 
equation (2.6) and separating real and imaginary parts give 

n t + div (n VS) = 0, (2.39) 

St + ||V5-| 2 + nn + \ (x> + 7 > 2 + 7z V) = t _L Av ^. (2.40) 

The equation (2.39) is the transport equation for the atom density and (2.40) the 
Hamilton- Jacobi equation for the phase. 

By formally passing to the limit e — > (cf. [30]), we obtain the system 

n° t + div (n° VS°) = 0, (2.41) 
S o + 1 )V5 0|2 + Kn o + 1 ^2 + 7 y + ^ = Q {2A2) 

It is well known that this limit process is only correct in the defocusing case k > 
before caustic onset, i.e. in time-intervals where the solution of the Hamilton- 
Jacobin equation (2.40) coupled with the atom-number conservation equation (2.39) 
is smooth. After the breakdown of regularity oscillations occur which make the 
term y A^/n at least 0(1) such that the validity of the formal limit process is 
destroyed. The limiting behavior after caustics appear is not understood yet except 
in the one-dimensional case without confinement, see [38]. Also, the focusing case 
k < is not fully understood yet. 



3 Numerical approximation 

In this section we present a time-splitting Fourier spectral method, which was used 
by Bao et al. to numerically solve the Schrodinger equation in the semiclassical 
regime [7, 8]. We reiterate that neither time splitting discretisations nor Fourier 
spectral methods are new, both have been applied successfully to many PDE prob- 
lems [13, 32, 54]. Here we adapt the combination of both techniques to the GP 
equation and infer computational domain and mesh size controls from analytical 
(perturbation) results. The merit of this method is that it is unconditionally stable, 
time reversible, time-transverse invariant, and conserves the total particle number. 
Also, it has very favorable properties with respect to efficiently choosing the spa- 
tial/temporal grid in dependence of the semiclassical parameter e. For simplicity of 
notation we shall introduce the method in one space dimension (d = 1). General- 
izations to d > 1 are straightforward for tensor product grids and the results remain 
valid without modifications. For d — 1, the equation (2.34) with periodic boundary 
conditions becomes 

ie — ^— = -—ip xx {x,t) + yVOM) + K i\^{x,t)\ 2 i)(x,t), a<x < b, (3.1) 

i{j(x, t = 0) = i/>°(x), a<x <b, (3.2) 
rj}{a,t)=tj}{b,t), Ma,t)=Mb,t), t > 0. (3.3) 
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We choose the spatial mesh size h = Ax > with h = {b — a)/M for M an even 
positive integer, the time step k = At > and let the grid points and the time step 
be 

Xj := a + j h, t n := n k, j = 0, 1, • • • , M, n — 0, 1, 2, • • • 

Let ip™ be the approximation of ip(xj,t n ) and ip n be the solution vector with com- 
ponents ■?/>". 

Time-splitting spectral method (TSSP). From time t — t n to t — t n+1 , the 
GPE (3.1) is solved in two splitting steps. One solves first 

e 2 

lei^t = -—ip xx , (3.4) 
for the time step of length k, followed by solving 

i £ ^iM = ^( x ,t) + Kl \^(x,t)\ 2 ij(x,t), (3.5) 

for the same time step. Equation (3.4) will be discretized in space by the Fourier 
spectral method and integrated in time exactly. For t e [t n ,t n+ i\, the ODE (3.5) 
leaves invariant in t [7, 8] and therefore becomes 

i £ ^jM = ^( x ,t) + KMx,t n )\ 2 tP(x,t) (3.6) 

and thus can be integrated exactly. From time t — t n to t — t n+1 , we combine the 
splitting steps via the standard Strang splitting: 



iP* = e -i(^ 2 /2+«i|^| 2 )fc/(2e) 
Af/2-1 



J 



^T = T7 E e-^V 2 $ e^~ a \ j = 0, 1, 2, • • • , M - 1, 

M l=-M/2 

= e -i^/2 +Kl \^)k/(2e) ^« j = 0, 1, 2, • • • , M — 1, (3.7) 

where the Fourier coefficients of are defined as 

2tt/ ~ ^ ... . M M 



Hi = 7^, ft = E # e-^- a \ / = -_..., y - l. (3.8) 



The overall time discretization error comes solely from the splitting, which is second 
order in k for fixed e > 0. The spatial discretization is of spectral (i.e. 'infinite') 
order of accuracy for e > fixed. An error analysis for linear Schrodinger equations 
taking into account the e-dependence of the global error can be found in [7] (where 
it is shown that k = 0(1) and h = ^ = 0(e) give correct observable), numerical 
tests for nonlinear problems in the semiclassical regimes in [8]. More restrictive 
meshing strategies are typically necessary in nonlinear cases, cf. Section 4.1. 
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For comparison purposes we review now alternative numerical methods [2, 16, 26] 
which are currently used for solving the Gross-Pitaevskii equation of BEC. One is 
the Crank-Nicolson finite difference (CNFD) scheme [2]: 



n+l 
3 



k 



I X 



ie 

2 



•n+l 



• n+l 



4e 



3 + 



C +1 = V'm , 



Ye 

J,n+l _ in+1 
YM+l — Wi , 

M. 



j = l,2,-..,M, 



'.n+l 

tf = Mxj), j = o,i, 

Another one is the Crank-Nicolson spectral (CNSP) method 
W +1 - M 



k 



IE 

7 



+i 



IT 



^n+l _ ^ 



^ = M^j), j = 0,l,-..,M; 
where D^ x , a spectral differential operator approximation of d xx , is defined as 



dLu 



M/2-1 
l=-M/2 



(3.9) 



Both methods are unconditionally stable, time reversible, conserve the total particle 
number but they are not time transverse-invariant. We do not present comparism 
tests with fully implicit and fully explicit finite difference methods since they are 
not at all competitive with the time splitting-spectral method. Generally 1) they 
require severe stability constraints on the mesh sizes, 2) they do not conserve the 
total particle number, 3) they are not time transverse invariant. For a mathematical 
analysis of FD-methods for Schrodinger type equations in semiclassical regimes we 
refer to [43, 44]. 



4 Numerical examples 

In this section, we first perform a numerical comparison of TSSP, CNFD and CNSP 
in terms of accuracy and mesh size strategy for a Id defocusing GPE. Then we 
apply the TSSP for solving Id, 2d and 3d GPEs of Bose-Einstein condensation. 
Furthermore we also give a physical discussion on our numerical results. 

In our computations, the initial condition for (2.34) is always chosen in WKB 
form: 

^(x, t = 0) = V°(x) = A°(x) e ,s °W/ £ , (4.1) 

with A and S° real valued, regular and with A°(x) decaying to zero sufficiently 
fast as |x| — > oo. We compute with TSSP on a domain, which is large enough (as 
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controlled by the initial data) such that the periodic boundary conditions do not 
introduce a significant aliasing error relative to the whole space problem. There are 
certainly more sophisticated analysis for controlling aliasing errors, however these 
do not significantly improve the results for exponentially decaying initial densities. 

To quantify the numerical results, we define the condensate widths along the x, 
y and z-axis as 

<r* = y/((x-(x))*), cr y = y/((y-(y))*), a z = J((z-(z)y), (4.2) 
where brackets denote space averaging with respect to the position density 

(/>= f /(x)|^(x,t)| 2 rfx. 

JR d 

4.1 Comparisons of different methods 

Example 1 Id defocusing condensate, i.e. we choose d — 1 in (2.34) with positive 
K\. The initial condition is taken as 

VKz, 0) = -^e-* 2 /( 2£ \ .61 (4.3) 

We choose e = 0.1 and K\ = 1.2649 and solve this problem on [—16, 16], i.e. a = — 16 
and b = 16 with periodic boundary conditions. Let ip be the 'exact' solution which 
is obtained numerically by using TSSP with a very fine mesh and time step, e.g., 
h — and k = 0.00001, and iph,k be the numerical solution obtained by using a 
method with mesh size h and time step k. 

First we compare the discretization error in space. We choose a very small 
time step, e.g., k = 0.00002 such that the error from the time discretization is 
negligible compared to the spatial discretization error, and solve the GPE using 
different methods and varying spatial mesh sizes h. Table 1 lists the numerical 
errors \\ip{t) — ^h,k(t)\\i 2 at £ = 2 for varying spatial mesh sizes h. Clearly TSSP and 
CNSP show roughly the same errors due to the fact that the temporal discretization 
is almost 'exact'. 

Secondly, we test the discretization error in time. Again, we take e — 0.1 and 
Ki = 1.2649. Table 2 shows the numerical errors \\ip(t) — iph,k{t)\\p at t = 2 with 
a very small mesh size h = ^ for different time steps k and different numerical 
methods. Here, CNSP and CNFD show almost no difference since the spatial dis- 
cretization now is almost 'exact'. 

We also tested numerically the unconditional stability of the time-splitting spec- 
tral method, which was already proven rigorously in [7] . Numerical tests showed no 
significant accumulation of round-off errors and conservation of the discrete Z 2 -norm 
was observed up to 10 significant digits for tests performed with h — ^ and k = 0.2, 
k = 0.05, k = 0.01 computing up to t = 4. 
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Mesh 




h = i 




h = ^ 




lb 32 




CPU time 


TSSP 


0.2248 


2M8E - 


2 


3.64LE- 


5 


7.982E - 


10 


0.01s 


CNSP 


0.2248 


2.048E - 


2 


3.642E - 


5 


8.538£ - 


8 


15.54s 


CNFD 


0.6314 


0.3380 




8.784E - 


2 


2.801£ - 


2 


1.27s 



Table 1: Spatial discretization error analysis: \\4>(t) — ^h,k{t)\\i 2 at time t = 2 
under k = 0.00002. The CPU time is counted at the same accuracy (i.e. ||-0(2) — 
^,fc(2) ||p 3.65E—5) and on an AlphaServer DS20 workstation. For that accuracy, 
TSSP needs h = ^ and k = 0.001, CNSP needs h = ^ and k = 0.00002 and CNFD 

Id Id 

needs h = -r^rj and k = 0.00001 



Time step 


k = 0.05 


k = 0.025 


k = 0.0125 


fc = 0.00625 


TSSP 


1.112E-2 


1.716S-3 


4.021E-4 


1.045E - 4 


CNSP 


0.5215 


0.1247 


4.363£ - 2 


1.565S-2 


CNFD 


0.5344 


0.13720 


6.12LE - 2 


3.723.B - 2 



Table 2: Time discretization error analysis: \\ip(t) — iphkif)\\i 2 at time t = 2 under 
h = ± 

n 32' 

At last, we test the e-resolution of different methods. Here we shall compare 
the meshing strategies required in order to get the 'correct' condensate density 
\ip\ 2 , for different methods when decreasing the semiclassical parameter e. Figure 
1 shows the numerical results with different combinations of e, h, k for different 
methods. Furthermore Figure 2 shows the evolution of p = \ip\ 2 in space-time 
and the condensate width as a function of time by using TSSP for e — 0.1 and 
«i = 1.2649. 

From Tables 1-3 and Figure 1, one can make the following observations: 

(1) . For TSSP, the spatial and temporal discretization errors are of spectral and 
second order accuracy, respectively. The admissible meshing strategy for obtaining 
the 'correct' condensate density in the defocusing case is: h = 0(e) and k = 0(e). 
This method is explicit, unconditional stable and its extension to 2d and 3d cases 
is straightforward without additional numerical difficulty 

(2) . For CNSP, the spatial and temporal discretization errors are also of spectral 
and first order accuracy, respectively. But the admissible meshing strategy is: h = 
0(e) and k = o(e). Furthermore this method is implicit and its extension to the 2d 
or 3d case is expensive except when an ADI technique is used, which destroys the 
spatial spectral accuracy. 

(3) . For CNFD, the spatial and temporal discretization errors are of second 
and first order, respectively, and the admissible meshing strategy is: h = o(e) and 
k = o(e) (see [43]). This method is implicit and the remark of (2) applies. 

Furthermore, the storage requirement of TSSP is less than the other two meth- 
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ods. The number of operations needed per time step is 0(M In M) for TSSP, at 
least 0(M 2 ) for CNSP, and O(M) for CNFD when an ADI technique is used in 2d 
and 3d, where M is the total number of unknowns. To attain the same order of 
accuracy, CNFD needs many more grid points than TSSP. 

4.2 Applications 

Example 2 2d defocusing condensate, i.e. we choose d = 2 in (2.34). We solve this 
problem on [— 8,8] 2 with mesh size h = ^ and time step k = 0.001. We present 
computations for four cases: 

/. 0(1) -interactions, zero initial phase data 

£ = l.o, ly = 1.0, k 2 = 2.0 (7 2 = 10.0, 5 = 1.586), 
#r,y,0) = _L e -(* 2 +s/ 2 )/(2e). 

\/7T£ 



II. Very weak interactions, anisotropic condensate, nonzero initial phase 
£ = l.o, ly = 2.0, k 2 = 0.1 (7 2 = 10.0 , 5 = 0.0793), 



7T£ 



Strong interactions, nonzero initial phase data 

£ = 0.1, 7y = 1.0, f£ = ^K 2 ^ y /n, k 2 = 1.2649 (j z = 10.0, 5 = 65.5227), 

^,y,0) = ( ^(^-(^ + y")/2)/ K2 e COsh ( V ^) /£ , x 2 + < 2^, 
[ 0, otherwise. 

IV. 0(1) -interactions, anisotropic condensate with changing trap frequency 
£ = 1.0, ly = 2.0, e x = 2.0, k 2 = 2.0 (7 2 = 10.0, 5 = 1.586), 

1/4 

tf,(x,y,0) = J^ e -(- 2 +7^ 2 )/(2 £ i). 

x /7T£i 



Figure 3 shows the surface plot of p = \ip\ 2 (labeled as \u\ 2 in the figures) at time 
t — 40 and the condensate widths as a function of time for case I. Furthermore, 
Figure 4 shows similar results for case II, Figure 5 for case III, Figures 6-7 for case 
IV. 

Example 3 2d focusing condensate, i.e. we choose d = 2 in in (2.34). We solve 
this problem on [—10, 10] 2 with mesh size h = ^ and time step k = 0.00005. We 
present computations for three cases: 
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7. 0(1) -interactions, positive initial energy 

£ = 1.0, ly = 1.0, k 2 = -2.0 (7 2 = 10.0, 5 = -1.586), 
^(x,y,0) = ^e-^ + y 2 ^ 2£ \ 

\/7T£ 



77. Strong interactions, negative initial energy 

e = 0.3, 7 y = 1.0, k 2 = -1.9718 = 10.0, 5 = -7.545), 



#r,y,0) 



1 



-(x 2 +w a )/(2e) 



Figure 8 shows the surface plot of p = \ip\ 2 at time t — 40 and the condensate 
widths as a function of time for case I. Furthermore, Figure 9 for case II. 

Example 4 3d defocusing condensate, i.e. we choose d = 3 in in (2.34). We 
solve for x e [—8, 8] 3 with mesh size h — | and time step k = 0.001 and present 
computations for two cases: 

I. Anisotropic condensate with changing trap frequency 

e = 1.0, % = 2.0, 7, = 4.0, e\ = \, «3 = 0.1 (8 = 0.1), 
0) = ( ^) V4 e -(^ +7y ^+T^ 2 )/(2 £l ) - 

77. Cylindrically symmetric condensate with changing trap frequency 

e=1.0, 7„ = 1.0, 7, = 2.0, £i = ^, k 3 = 1.0 (5 = 1.0), 

0) = ( 7" 72)1/4 e -(^+7^ 2 +7^ 2 )/(2ei)^ 

'(tci) 3/4 



Figure 10 shows the condensate widths as a function of time for cases I and II. 

Example 5 2-dim vortices in BEC. We choose d = 2 and simulate the effect 
of stirring the (stable) condensate by adding a narrow, circularly moving Gaussian 
potential W(x, t) to the stationary trap potential in (2.34). W(x,t) represents, for 
example, a far-blue-detuned laser [14]. We set 



W(x,t) = W s (t)exp -4|x-x s (t)| 2 /V fl 



with the center x s (£) = (r cos u s t,r sin u s t) moving on a circle with radius r and 
frequency u s . We start the simulation with the ground state in 2 dimensions (no 
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stirrer at t — 0) and minimize transient effects by increasing the stirrer amplitude 
W s (t) linearly from at t = to a final value W s (t = ir) =: Wf at t = ir. The 
stirrer is then linearly withdrawn from t = 4n to t = 5n (after constant stirring, 
i.e. W s (t) = Wf for 7r < t < An) and the condensate is left to evolve freely after 
t = 57r. We recall that 2-dim and 3-dim vortices simulations were already performed 
in [14, 15, 3], here we present this example in order to put our numerical method to 
an important physical test. 

We take the numerical values e = l/y/50, n d = 1, ^ y = 1, Wf = y/2, V s = 1/50 1 / 4 , 
ro = 2/50 1//4 and uo s = 1 for our simulation. We remark that a vortex in the 
condensate is a point x„ with i/j(x. w ) = and singular or undefined phase. 

In Figure 11, we show the contour plot of the density |t/>(x, t)| 2 at t — Ylix (where 
the fluid has already settled down after the stirring) and x, ^-sectional plots of the 
vortex centered at (x — 0.141, y w —0.229). in fact three vortices (labeled by 'X'), 
located at (-0.141, -0.229), (1.093, -0.0353) and (0.282, 1.481) were identified. For 
an analysis of vortex-formation in semiclassical limits of the Schrodinger equation 
we refer to [42]. 

4.3 Discussion 

In Sec. 4.1 we compared different numerical methods for solving the GPE with the 
TSSP. Now we complete our investigation on the validity of using the TSSP for 
solving the GPE by comparing the results obtained in Sec. 4.2 with well known 
properties of Bose- Einstein condensates at very low temperatures. 

In example 1 we present Id simulations. Initially the condensate is assumed to 
be in its noninteracting ground state when at t — repulsive interaction is turned 
on. In current experiments a change in the interaction strength can be achieved by 
applying external magnetic fields. Close to a Feshbach resonance the interaction 
strength shows a strong dependence on the magnetic field and even its sign can be 
changed by an appropriate choice of the magnetic field [49]. At the same time we 
also change the trap potential by setting e = 0.1. These sudden changes lead to 
many rapid oscillations in the condensates (cf. Fig. 2a)). The dominant excitation 
caused by the interactions is, as expected [37, 25], the oscillation of the condensate 
width at approximately twice the trap frequency (cf. Fig. 2b)). 

Example 2 presents 2d simulations for various cases. In cases I there are 0(1) 
- interactions and in case II we investigate a weakly interacting condensate. We 
assume the condensate to be initially in the noninteracting ground state with pos- 
sibly a nonuniform phase. As in example 1 turning on the interactions causes the 
condensate to oscillate at approximately twice the corresponding trap frequency. 
Higher excitations like in Fig. 2 a) are not visible in Figs. 3 a), 4 a) since compared 
to example 1 the strength of the interactions is much smaller. A nonuniform phase 
like in II causes the amplitude of the oscillations to be different in x and y direc- 
tion. Note also that while in case I the condensate immediately starts to expand 
(cf. Fig. 3 b)) (due to the repulsive interactions) it starts to contract for the initial 
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condition with a nonzero phase in case II. In case III we investigate the evolution 
of a strongly interacting condensate initially in the (approximate) Thomas Fermi 
ground state with an additional phase. Again, since this is not the exact ground 
state (see Remark 2.2) the width of the condensate starts to oscillate at about twice 
the trap frequency (cf. Fig. 5a)). In this case, however, due to the strong nonlin- 
earity the oscillations along the x and y direction are coupled with each other as 
can be seen from Fig. 5 b). In case IV we investigate the effect of changing the trap 
frequency and turning on repulsive interactions (see Fig. 6). Like in the previous 
cases we find the dominant effect to be oscillations at about twice the corresponding 
trap frequency. The condensate initially starts to contract since the trap frequencies 
are increased at t — 0. For the initial conditions chosen in this case the amplitudes 
are sufficiently large to swap the widths a x and a y (cf. Fig. 7), whereas for small 
changes in the trap frequency no swapping of the condensate widths would occur. 
The numerical results for the oscillations of a BEC obtained in example 2 agree very 
well with experimental and theoretical results [37, 25]. 

In example 3 we show solutions for a focusing nonlinearity in 2d. Case I shows 
the effect of turning on 0(1) attractive interactions which leads to oscillations similar 
to those discussed in the previous examples (see Fig. 8). In case II a condensate 
with negative initial energy is shown. We have not discussed the reduction of the 
GPE to 2d for this case. Also, our simulations do not contain loss terms which 
become important in condensates at high densities. Thus we do not give a physical 
interpretation of these results. However, this example shows that the numerical 
method is applicable to the case of strong focusing nonlinearities in the GPE. Our 
numerical results confirm that the attractive GPE in 2d with negative initial energy 
will blow up at finite time (cf. Fig. 9). Furthermore, we point out that the TSSP 
allows the inclusion of loss terms into the GPE and it is also feasible to solve the 
GPE in 3d for the attractive case [10]. Therefore the TSSP is a promising candidate 
for simulating the recent experiments on collapsing and exploding BEC's by Donley 
et al. [24] which requires full 3d simulations and the inclusion of loss channels. 

Example 4 shows the effects of turning on repulsive interactions and changing 
the trap frequency in a 3d condensate. As expected from our previous simulations 
we see in Fig. 10 oscillations at twice the trap frequency in directions x, y and z, 
respectively. The amplitude of the oscillations decreases with increasing frequency, 
i.e. it becomes more difficult to excite oscillations for larger trap frequencies. This 
behavior is one of the basic assumptions allowing the reduction of the GPE to 2d 
and Id in the cases where one or two of the trap frequencies are much larger than 
the others (cf. Sec. 2.3). 

The last example 5 shows the creation of vortices in a 2d BEC by stirring it with 
a blue detuned laser beam (see also [14, 15, 45]) where we identify the creation of 
three vortices in the BEC as shown in Fig. 11. We note that since we study the GPE 
without taking into account the interaction between the condensate and a thermal 
cloud of atoms (additional dissipation) no stationary state showing an Abrikosov 
lattice of vortices is found as in recent experiments [45] and numerical studies on 
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the effects of a thermal cloud [47] on the vortex formation. However, we point out 
that full 3d simulations including the effects of thermal particles with a very high 
precision are feasible based on the TSSP. 

Finally, we note that the TSSP is a very powerful versatile numerical method 
for solving the GPE which can be applied to a large number of different physical 
situations. The efficiency of this method and the high precision of the solutions make 
the TSSP a good choice for solving experimental situations that are numerically very 
demanding. Among these we believe that numerical studies of collapsing condensates 
with attractive interactions and multi-component condensates taking into account 
all experimentally relevant extensions of the GPE as well as on extensions of the 
GPE dealing with dissipation mechanisms [31] will be feasible by using the TSSP 
[10]. 

5 Summary 

We studied a numerical method for solving the time-dependent GPE which de- 
scribes trapped Bose-Einstein condensates at temperatures T much smaller than 
the critical condensate temperature T c . We started with the 3d GPE, scaled it to 
obtain a four-parameter model, and showed how to approximately reduce it to a 2d 
GPE and a Id GPE in certain limits. We provided the approximate ground state 
solution of the GPE in two extreme regimes: (very) weakly interacting conden- 
sates and condensates with strong repulsive interactions. Then, most importantly, 
we used the time-splitting spectral method in connection with analytical considera- 
tions based on perturbation theory (mesh-size control, dimension reduction) to solve 
the time-dependent GPE in Id, 2d and 3d. Extensive numerical examples in Id, 2d 
and 3d for weakly /strongly interacting condensates defocusing/focusing nonlinear- 
ity, and zero/nonzero initial phase data were presented to demonstrate the power 
of the time-splitting spectral numerical method. Finally, we want to point out that 
equations very similar to the GPE are also encountered in nonlinear optics. In the 
future we plan to apply this powerful numerical method to physically more complex 
systems like multi component condensates, collapsing condensates with attractive 
interactions and also to describe coherent atomic samples in wave guides. 
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Figure 1: e- resolution comparison in Example 1 for condensate density \ip\ 2 of differ- 
ent methods, ' — ': 'exact' solution, '+ + +': numerical solution, a). TSSP: e = 0.4, 
h=\,k = 0.04 (left); e = 0.1, h = ±, k = 0.01 (middle); and e = 0.025, h = ^, 
k = 0.0025 (right), b). CNSP: e = 0.4, h = \, k = 0.02 (left); e = 0.1, h = ±, 
k = 0.005 (middle); and e = 0.025, h = ±, k = 0.00125 (right), c). CNFD: e = 0.4, 
h = jq, k = 0.02 (left); e = 0.1, h = ±, k = 0.005 (middle); and e = 0.025, h = ±, 
k = 0.00125 (right). 
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Figure 2: Numerical results in Example 1 for e = 0.1 and Ki = 1.2649. a). 
Evolution of the position density; b). widths of the condensate as a function of 
time. 
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Figure 3: Numerical results in Example 3 for case I. a). Surface plot of the position 
density at t — 40; b). widths of the condensate as a function of time. 
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Figure 4: Numerical results in Example 3 for case II. a). Surface plot of the 

position density at t — 40; b). widths of the condensate as a function of time. 
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Figure 5: Numerical results in Example 3 for case III. a). Surface plot of the 

position density at t — 40; b). widths of the condensate as a function of time. 
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Figure 6: Numerical results in Example 3 for case IV. a). Surface plot of the 

position density at t — 6.8; b). widths of the condensate as a function of time. 
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Figure 8: Numerical results in Example 4 for case I. a). Surface plot of the position 

density at t — 40; b). widths of the condensate as a function of time. 
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Figure 9: Numerical results in Example 4 for case II. a). Surface plot of the 
position density at t — 0.5; b). peak of the position density \ip(0, 0, t)\ 2 as a function 
of time. 
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Figure 10: Widths of the condensate as a function of time in Example 5. a). For 

case I; b). for case II. 
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Figure 11: Numerical results at time t = Yin in Example 6. a). Contour plot of 
the density function \ip\ 2 (three vortices were identified and labelled by 'X'); b). x, y- 
sectional plot at a vortex (center of it is labelled by 'O') located at (—0.141, —0.229); 
'— ': \i/j(x, -0.229, 12tt)| 2 , '- - - ': |^(-0.141, y + 0.088, 12vr)| 2 . 
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